Estimation of Initial Conditions and Secure Communication 
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We estimate the initial conditions of a multivariable dynamical system from a scalar signal, 
using a modified Newton-Raphson method incorporating the time evolution. We can estimate initial 
conditions of periodic and chaotic systems and the required length of scalar signal is very small. We 
also find that the information flow from one variable to the other has logarithmic dependence on 
time. An important application of the method is in secure communications. The communication 
procedure has several advantages as compared to others using dynamical systems. 
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A trajectory of a given dynamical system in its state 
space depends on the set of initial conditions with which 
it starts. In particular, the state of a chaotic system at 
a later time is exponentially sensitive to changes in its 
initial state ||] . This feature of a chaotic system leads to 
a complex behaviour in state space that appears random 
yet is deterministic, and the time evolution is uniquely 
fixed by the initial state of the system. Though there are 
several invariant measures of a chaotic system which are 
not sensitive to the initial conditions, the exact trajectory 
crucially depends on the initial state and hence is difficult 
to reproduce due to sensitivity to initial conditions. 

In light of these facts, it is interesting and important 
to ask whether the complete set of initial conditions of a 
given multivariable dynamical system can be estimated 
from a given scalar time series of a single state space vari- 
able. In this letter we present a novel and simple method 
to estimate the initial conditions from a given scalar time 
series. The method is based on a modified multidimen- 
sional Newton-Raphson method that includes the 
time evolution of the system. The length of the time se- 
ries required for estimating the initial conditions is very 
small. Also, the method works even when the conditional 
Lyapunov exponents are positive. 

An important application of our method is in the area 
of secure communications A nice feature of this 

application is that the signal that is transmitted is not 
the one modulated by the information signal, making it 
difficult to crack the method. Our method is also useful 
in the problem of synchronization of chaotic signals 
The knowledge of initial conditions means the response 
system can be be synchronized with the drive system al- 
most instantaneously thereby removing the problem of 
transients. Also synchronization can be achieved in most 
of the cases where other methods fail ||^. 



Let us consider an autonomous dynamical system 
given by, 



F(x), 



(1) 



where x = (xi,X2, ■ ■ ■ ,Xd) is a d-dimensional state vec- 
tor whose evolution is governed by the function F = 
{Fi,F2, ■ ■ ■ ,Fd)- Given an initial state vector x(0) at 
time i = 0, the time evolution x(t) is uniquely deter- 
mined by Eq. . Now let us assume that only one com- 
ponent of the state vector is known to us and we take it 
to be xi{t) without loss of generality. The problem that 
we address is to obtain the initial state vector x(0) from 
the knowledge of the scalar signal xi{t). 

Let y(0) denote a random initial state vector and y{t) 
its time evolution obtained from Eq. (|l|). Let w(t) de- 
note the difference w(t) = y{t) — x(t). We look for the 
solution of the equation 



w(t) = 0. 



(2) 



Noting that the initial state vectors y(0) and x(0) 
uniquely determine the difference w(i), one of the solu- 
tions of Eq. (H) is y(0) — x(0) = and this is the solution 
that we are searching for. 

We now introduce the notation w" = w"(y°,x") = 
w(nAi), where At is a small time interval. Similarly, 
y" = y{nlS.t) and x" ~ x(nAt). With this notation 
condition (||) can be written as w" — 0. 

Our approach to the solution of Eq. (||) is a modified 
Newton-Raphson method which includes the time evo- 
lution of the system. 

Let us first consider w^. We have 

= wi(x°,x°) = wi(y"+<5y°,x°), 
= wi(yO,x") + (5yO.Vyo)wi(yO,x°) + 0((<5y°)2), (3) 
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where (5y° = — y° = — w" and the last step is a Taylor 
series expansion in Sy'^. For small Ai, we can write 

wi(y°, x") = w° + A< [F(y") - F(x°)] + ©((At)^). (4) 

Substituting Eq. in Eq. (^) and neglecting higher or- 
der terms, we get, 

wi(y°,x°) = w° - Ai(w°.Vyo)F(y°). (5) 

It is convenient to write the above equation in a matrix 
form as 

= (I + At J°)W° = A°W° , (6) 

where W"^ is the column matrix corresponding to the 
vector w^, I is the identity matrix, A" = / + At J", 
and the elements of the Jacobian matrix J" are J," = 

Next we consider or W^. Proceeding as above, we 
get (see Eq. (|)), 

W^ = {I + At J^)W^ = (I + At J^){I + At J°) W° 
= A^A°W°. (7) 

Similarly, the equation for is 

M/" = A"-M"-2...^0p^o^ (8) 

We now concentrate on the first component of the sig- 
nal whose time series is assumed to be known. For a 
d-dimensional system we need d — 1 equations to deter- 
mine the initial state vector x°. Eqs. (^), (|^) and (||) give 
us the required relations. 

i=i 

d 

wt'= E 4r'---^LM^^, (9) 

i,...,^,m— 1 

These are d~ \ simultaneous equations for W'^ . 

The numerical procedure is as follows. We set the ini- 
tial state of system (|l|) to a random initial guess vector 
(y°)o!d ^'^^^ iyi)oid = ^1 evolve it using Eq. (|). 
Using this vector y{t) we write down d — 1 simultaneous 
equations (Eqs. (^)) which can be solved for d — 1 un- 
known components of = — (5y°. Also, Syl — 0. Thus 
the initial guess vector can be improved by 

(y°)„e. = (y°)o.. + ^y°- (lo) 

This sets up an iterative scheme giving us better and 
better estimates of the initial vector which converge to 



We note that as in Newton-Raphson method, the 
choice of the initial guess vector can be very impor- 
tant 1^. In some cases, the iterative procedure of Eq. ( p^ ) 
may not converge or converge to a wrong root. In such 
cases, a different choice of initial guess vector can be use- 
ful. 

We further note the similarity of our method with the 
so called method of variational equations in analytical 
dynamics The method of variational equations can 
be applied to a known Hamiltonian system to determine 
an unknown neighboring trajectory to an already known 
one. There, the method requires a complete particular 
solution of a known set of Hamiltonian equations of mo- 
tion. In contrast, we have used our method for dissipative 
chaotic systems. In such systems an analytical solution 
of the equations of motion cannot be known. Further 
our method requires only one component of a complete 
trajectory to be sampled. This is important for the ap- 
plication of our method to secure communications as we 
will demonstrate afterwards. 

We now illustrate our method of estimating the initial 
state. As an example we discuss the Rossler system given 
by 0, X = (-X2 - X3,xi + ax2jb + X3,[xi - c)). First, 
we consider a case when the time series for xi is given 
and we want to estimate (a;^,^^). We chose the parame- 
ters (a, 6, c) such that the system is in the chaotic regime 
and the initial state x° is on the chaotic attractor. We 
start with an arbitrary initial state y" — {ui, 1/2, y 3) with 
yi — Xi- From Eqs. (^ we get a pair of simultaneous 
equations for {6y2,Sy^), which can be solved to obtain 
these values. With 5y'^ = we use these in an iterative 
manner (Eq. (p^)) to obtain the correct initial conditions. 

Let Ci = \y^ — x'^l denote the absolute error in the es- 
timation of a:°. In Fig. 1(a) we plot a graph of errors 62 
and 63 plotted against the number of iterations, n, of our 
method (Eq. (pO|)). From Fig. 1(a) we see that the errors 
go to zero and the successive estimates converge to the 
correct values of (a;2,a;3). Using only two data points in 
the given time series Xi{t), we can thus readily estimate 
the full initial state x°. We also note that the rate of 
convergence is very good. In about 8 to 10 iterates we 
obtain the initial values (a;^, x'^) to within computer accu- 
racy. If we write the deviations of the successive iterates 
from the correct values in the form 

(e.)„ = |(y°)„-^°he-"". (11) 

where n is the number of iterations, then the value of 
the parameter a is found to be 2.12 for 62 and 2.10 for 
63 . This is consistent with the fact that Newton-Raphson 
method has a quadratic convergence ||l|,||. 

We note that the largest Lyapunov exponent for the 
subsystem (y^iJ/s) (conditional or subsystem Lyapunov 
exponent) is positive [pi . The success of our method does 
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not depend on whether this Lyapunov exponent is pos- 
itive or negative. This is important for synchronization 
of chaotic signals. 

We next present cases where time series for the vari- 
ables x-2 and of the Rossler system are given. The 
procedure is similar to the case of time series for as 
discussed above. Fig. 1 (b) shows the errors ei and 63, 
when time series for X2 is given, plotted against the num- 
ber of iterations. The parameter a (Eq. ([TTI)) is 2.08 for 
ei and 2.10 for 63. This again indicates a quadratic con- 
vergence. Similarly, Fig. 1 (c) shows the quantities ei 
and 62 when time series for x^ is given, as a function of 
the number of iterations. The parameter a (Eq. jTl])) is 
2.02 for d and 2.10 for 62- We note that the largest sub- 
system Lyapunov exponent is negative when time series 
for X2 is given and is positive when time series for X3 is 
given H. 

Wc have successfully applied our method to estimate 
the initial state vector using a given scalar time series for 
many other dynamical systems as well. These include 
Lorenz system in its periodic, chaotic or intermittent 
regimes, Chua's circuit j^], the disk dynamo system mod- 
eling a periodic reversal of earth's magnetic field JlO| , pl| , 
a 3-d plasma system formed by a three wave resonant 
coupling equations [T^ and a four dimensional phase con- 
verter circuit psf . 

We now turn to the important question of the rate of 
flow of information about other variables into the given 
variable, i.e. the variable whose time series is known, say 
xi. To investigate this question we start with the initial 
value a:i(0) at t = 0. This alone does not tell us anything 
about the other variables. As time evolves information 
about other variables flows into Xi and the method pre- 
sented above allows us to extract this information. As a 
specific example consider the Rossler system with time 
series for Xi . Fig. 2 shows the plot of the asymptotic error 
62 of the estimation of 0:2(0) for large n, denoted by 600, 
as a function of the total time interval r = (d — l)At for 
which the time series of xi is used in the calculations. We 
observe that the accuracy of estimation improves with r 
and shows a power law behaviour of the type 



(12) 



Fig. 2 shows a very good straight line behaviour on a log- 
log plot over seven orders of magnitude and the slope 
gives the value of the exponent v = 2.00 ± 0.01. The 
power law dependence of the accuracy on time interval as 
in Eq. (12) shows that the information entropy in terms of 
bits or partitions grows logarithmically with time. Note 
that this flow of information is different from the usual 
linear growth of information entropy with time in chaotic 
systems which corresponds to an exponential dependence 
of 600 on T. The excellent straight line fit over seven or- 
ders of magnitude that we observe in Fig. 2, makes us be- 
lieve that the power law dependence of Eq. (O) and the 



corresponding logarithmic growth of information content 
with time is a genuine property of dynamical systems. 
Further, the exponent v appears to be universal and has 
the same value within numerical errors for the various 
examples that we have considered. 

Now we will discuss an important application of our 
method for secure communications. Let us suppose that 
Eq. (|l]) describes a c?-dimensional chaotic system at the 
transmitter end and a replica of the same at the receiver 
end. Let s(i), i = 1, 2, . . . , be the information signal to 
be encoded for communication and without loss of gener- 
ality let xi{t) be the scalar output signal to be transmit- 
ted. The signal s{t) is encoded as follows. While choosing 
the initial conditions of the vector x(0), the component 
xi(0) is chosen arbitrarily while the other d — 1 compo- 
nents are chosen as Xi(0) = s(i — 1), i = 2, . . . , d. Using 
the initial vector, x(0), we evolve Eq. (|]) for d — \ time 
steps of dt to obtain the first d values of the scalar vari- 
able xi. This is the first stage of evolution. The vector 
x(ti), where i„ = n{d— l)dt, is modified using s(i) as ex- 
plained below (Eq. (^3|)) and the resultant is evolved for 
d — 1 time steps completing the second stage. The pro- 
cess is repeated as often as required. After completion of 
each stage the vector x(t„) is modified as follows. 



Xi{tn) = Xi{tn) 

Xi{tn) = Xi{tn) + s{n{d - 1) + i - 1), i = 2, 



,d. (13) 



Thus, assuming that k, the number of elements of infor- 
mation signal s(z), is divisible by d — 1 we get k + 1 values 
of the encoded signal xi. 

The recovery of the information signal at the receiver 
end using our procedure of initial condition estimation is 
straightforward. At each stage d — 1 values of xi unable 
us to estimate the initial conditions thereby recovering 
the signal from Eq. ([l3|). 

The communication procedure described above has 
several advantages. Note that (a) The transmitted signal 
xi is not directly modulated by the information signal s. 
(b) The information signal is added to more than one 
variables. Due to these reasons, local approximation for 
the flow of the dynamical system in the embedded phase 
space [Q cannot be used to attack the coding proce- 
dure. Secondly, the transmitted signal is of the same size 
as that of the information signal, and there is no ini- 
tial transient as required in other synchronization based 
methods Thirdly, as in other methods using chaotic 
encoding the same signal s will give different outputs for 
reruns of the coding procedure thereby making any at- 
tack very difficult. 

Another important application of our initial condition 
estimation method is in the problem of synchronization 
of two identical chaotic systems coupled unidircctionally 
by a scalar signal. The knowledge of initial conditions 
enables us to obtain immediate synchronization thereby 
eliminating transients Also, the method could be 
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used repeatedly to maintain synchronization. This will 
be particularly useful for systems where the largest condi- 
tional Lyapunov exponent is positive and other methods 
of synchronization fail jsj . 

To summarize, we have introduced a novel yet sim- 
ple method to estimate initial conditions of a multivari- 
able dynamical system from a given scalar signal. Our 
method is based on a multidimensional Newton-Raphson 
method where we include the time evolution of the sys- 
tem. The method gives a reasonably fast convergence 
to the correct initial state. The required length of the 
time series is very small. The method works even when 
the largest conditional Lyapunov exponent is positive. 
We find that the information flow between different vari- 
ables grows logarithmically with time. This is different 
from the usual linear growth of information in chaotic 
systems. An important application of the method is in 
secure communications. This method of communication 
has several advantages as compared to other methods of 
communication based on nonlinear dynamical systems. 
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FIG. 1. Plot (a) shows the errors 62 and 63 on a logarithmic scale, as a function of n, the number of iterations of our method 
for Rossler system when time series for variable xi is given. The pluses show the values for e2 and the crosses those for 63. 
Errors are seen to approach zero as n increases. The parameters are {a,b,c) = (0.2,0.2,9.0), the time step At = 0.01 and we 
use fourth order Runge-Kutta method for the time evolution of Rossler equations. (We have also checked our results with fifth 
and sixth order Runge-Kutta methods.) Similarly plot (b) shows the errors ei (pluses), and eg (crosses) as a function of n 
when time series for X2 is given. Plot (c) shows the errors e| (pluses), and 62 (crosses) as a function of n when time series for 
a; 3 is given. 



FIG. 2. The points in the graph show the asymptotic error Coo of initial state estimation as a function of At. The case 
under consideration is the same as in Fig. 1(a) for 62 averaged over several realizations. A good straight line fit is obtained 
over seven orders of magnitude and has a slope of —2.0 ± 0.01. This indicates a logarithmic fiow of information from variable 
X2 to variable xi. 
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